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Abstract 

As a result of specific adaptations and habitat preferences strongly rheophilic fish species may show high levels of 
endemism. Many temperate rheophilic fish species were subjected to a series of range contractions during the 
Pleistocene, and then successfully expanded during the Holocene, colonising previously abandoned areas. The 
Carpathian barbel (Barbus carpathicus Kotlik, Tsigenopoulos, Rab et Berrebi 2002) occurs in the montane streams in 
three basins of the main Central European rivers in the northern part of the Carpathian range. We used genetic 
variation within 3 mitochondrial and 9 microsatellite loci to determine a pattern of postglacial expansion in 6. 
carpathicus. We found that overall genetic variation within the species is relatively low. Estimate of time to the most 
recent common ancestor (tMRCA) of mitochondrial sequences falls within the Holocene. The highest levels of 
genetic variation found in upper reaches of the Tisa river in the Danube basin suggest that glacial refugia were 
located in the south-eastern part of the species range. Our data suggest that the species crossed different 
watersheds at least six times as three genetically distinct groups (probably established in different expansion 
episodes) were found in northern part of the species range. Clines of genetic variation were observed in both the 
Danube and Vistula basins, which probably resulted from subsequent bottlenecks while colonizing successive 
habitats (south eastern populations) or due to the admixture of genetically diverse individuals to a previously uniform 
population (Vistula basin). Therefore, B. carpathicus underwent both demographic breakdowns and expansions 
during the Holocene, showing its distribution and demography are sensitive to environmental change. Our findings 
are important in the light of the current human-induced habitats alterations. 
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Introduction 

During the Pleistocene numerous temperate species were 
subjected to latitudinal range shifts and demographic 
fluctuations [Ij.The cooler climate and shortened vegetation 
period during glaciations radically altered the functioning of 
ecosystems. Species that were not able to establish 
populations in relic patches of the appropriate environment 
went extinct, while others that survived the colder period in 
glacial refugia recolonized temperate zones when the climate 
warmed. These range changes were determined by the ability 
of species to tolerate or adapt to changing conditions, and 
different species may thus have responded differently [2], [3j. 
In many species such retreat/expansion cycles were 
accompanied by demographic fluctuations [4j. 



Species dispersal depends strongly on the interconnectivity 
among patches of suitable environment and species with 
limited dispersal abilities thus may have been trapped in 
narrow patches of a suitable habitat. In riverine species, 
expansion corridors are limited to linear aquatic habitats within 
a river system. These pathways of expansion are isolated by 
watersheds, while in strongly rheophilic fishes migration may 
also be hindered by long stretches of lowland river segments 
[5]. 

The Carpathian barbel, Barbus carpathicus Kotlik, 
Tsigenopoulos, Rab et Berrebi, 2002, was only recently 
described as a separate species, and knowledge of life history 
specific to this species is therefore scarce [5], [6], [7]. The 
species occurs in Northern and Central part of Carpathian 
Mountains. Although some local populations of B. carpathicus 
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Table 1. Characterization of polimorphic microsatellite loci successfully amplified in the Carpathian barbel. 



Primer set 


Fluorescent dye 


Annealing temp. 


Locus 


Product size (bp) 


N alleles 


He 


Barb 37 


FAM-6 


56 


Barb 37 


201-224 


10 


0.43 


Barb 59 


TAMRA 


56 


Barb 59 


134-258 


41 


0.87 


Barb 79 


FAM-6 


54 


Barb 79 


200-216 


5 


0.53 


LC 293* 


Hex 


56 


LC 293b 


100-110 


5 


0.50 


LCo 4* 


Hex 


56 


LCo 4a 


207-211 


5 


0.13 








LCo 4b 


212-214 


3 


0.47 


MFW 17* 


FAM-6 


54 


MFW 17a 


180-190 


3 


0.14 








MFW 17b 


199-221 


9 


0.50 


MFW 19 


TAMRA 


54 


MFW 19 


187-297 


22 


0.80 



primers amplifying duplicated loci, 
doi: 10.1371/journal.pone.0082464.t001 



are large they are confined to montane/submontane streams 
and rivers [8]. 

The Carpathian barbel diverged from its closest congener 
Balkan barbel, B. balcanicus Kotlik, Tsigenopoulos, Rab et 
Berrebi, 2002, in the mid-Pliocene between approximately 3.9- 
4.1 million years ago (mya) [6]. Both species share many 
biological similarities with each other and with the third 
rheophilic barbel occurring in the Danube basin, Peteny's 
barbel, B. petenyi Kotlik, Tsigenopoulos, Rab et Berrebi, 2002. 
The similarities are such that the three above mentioned 
species have been confused as a single species until recently 
[7], yet there is a striking difference in the level of genetic 
variation between B. carpathicus and the other two species [6]. 
B. carpathicus shows great potential as a model for analysing 
the natural expansion of montane fish species. It belongs to a 
number of southern European species of montane barbels 
which share some common features. They are small to 
medium-sized, strictly rheophilic and therefore confined to 
montane ranges, and almost allopatric [9]. 

A recent study carried out in a small part of this area at the 
Vistula-Dniester watershed [5] based on variation within the 
nuclear microsatellite loci found evidence that (1) the cline of 
microsatellite diversity decreases eastwards, (2) the overall low 
genetic variation reflects probable population bottlenecks, and 
(3) that there are genetic similarities between populations in the 
adjacent rivers in the basins of the Vistula and Dniester. Such 
features of the populations of B. carpathicus north of the 
Carpathians resemble the gradual loss of genetic diversity 
which accompanies many postglacial expansions of species 
limited by hard barriers that permit only rare migrations [10], 
[11], [12], [13]. 

The aim of this study was to determine the pattern of 
expansion of B. carpathicus, i.e. (1) where the Pleistocene 
refugia were located, (2) whether B. carpathicus invaded rivers 
of the Baltic catchment's northern slopes of the Carpathians on 
a single or multiple occasions, (3) where the most likely points 
are at which the main Carpathian watershed was crossed, 
including whether the founders originated from a single or 
multiple sources, (4) what the likely origin of dines of genetic 
variation is, and (5) at what time these expansion events 
occurred, in particular whether the Vistulan populations were 
established in the Holocene or contain traces of earlier 



expansions or presence of glacial refugia. We addressed 
above question using sequence variation within mitochondrial 
haplotypes and length variation in nine nuclear microsatellite 
loci. 

Materials and Methods 

Ethics statement 

Standard procedures used in aquaculture for catching and 
marking fish were employed in the course of the study. No in 
vivo experiments were performed on animals. Permissions for 
catching and fin clipping were granted by the Voivode of the 
Matopolska region, Poland, Nos. RG. 1.6052-1/09, RG.I. 
6052-2/09, and RG.I.6052-6/09, Voivode of the Podkarpacie 
Province, Poland, Nos. RS.I.EK.6251/03/09 and RS.I. 
6251/08/08, Voivode of the Silesia Province, Poland, No: 
3905/TW/2009, Ministry of Environment of the Slovak Republic: 
Nos. 47/2008, 48/2008, 18/2009, 19/2009, and 6057/2010-2.1. 

No animals were sacrificed in the course of this study. The 
muscle samples were from previous published studies [6], [7] 
and were collected in 1997-1998. The sample collection was 
overseen by the Ethics Committee of the Institute of Animal 
Physiology and Genetics, Academy of Sciences of the Czech 
Republic. Fish were sacrificed using an overdose of 
anesthetics and muscle samples were taken from dead fish. 

Sample collection 

A total of 841 samples from 27 populations were analyzed. 
Sampling covered rivers from all three basins within the 
species range (Table 1, Figure 1), i.e. the Danube, Vistula and 
Dniester. To avoid analyzing single flocks, samples from 
several locations along the river were collected. Samples were 
pooled into distinct populations when sampling locations were 
divided by long segments of lowland river or the distance 
between sampling locations exceeded 100 km (e.g. Poprad vs. 
Dunajec). A unique sample collected in the Sikenica river 
belonging to the Hron basin was included in the Hron sample. 

The majority of samples (625) were fin clips collected in the 
years 2008-2009. Fish were immobilized with a portable 
pulsed-DC electroshocker, and after clipping 5-10 mm 2 of the 
tip of the pelvic fin they were released in the place of capture. 
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H10 H11 




Figure 1. Statistical parsimony network of mitochondrial haplotypes, and the distribution of haplotypes in current range 
of Barbus carpathicus. 

doi: 10.1371/journal.pone.0082464.g001 



Samples were preserved in 96% ethanol and stored in -20 °C 
upon arrival to the laboratory. 

We also used 66 fin clips collected for a previous study on 
the expansion of B. carpathicus [5]. To obtain muscle samples 
fish were euthanized using an overdose of 2-phenoxyethanol 
(bath in 0.4 mg/liter solution) without any prior manipulation. 

Laboratory protocols 

Immediately prior to DNA isolation, samples were desiccated 
overnight at 60 °C. Genomic DNA was isolated with 
NucleoSpin Tissue Kit (Macherey & Nagel), following the 
standard manufacturer's protocol. 

Three fragments of mitochondrial genome were amplified 
and sequenced in 78 individuals of B. carpathicus including the 
cytochrome b (cytB, 1193 bp): L14724 and H 1591 5 [14], NADH 
dehydrogenase subunit 2 (ND2, 1142 bp): ND2-2F (5'- 
CTTCCTTTACACCACTTTCT-3') and ND2-1R 

(5TTGAAGGCTTTTGGTCTAAT-3'), ATPase subunits 6 and 8 
(ATP, 809 bp): ATP8.2 and Co3.2 [15]. Both PCR conditions 
and thermal profiles were the same for the three loci. PCR 
reactions were prepared in 20 pi volumes containing 20 ng of 



template DNA, 200 nM of each primer (Genomed), 2 mM 
MgCI 2 , 200 nM dNTP, 1 x PCR buffer containing NH 4 S0 4 and 

1 U Taq DNA polymerase (Fermentas). Loci were amplified in 
PTC-200 MJResearch Thermal Cycler under the following 
thermal profile: 3 min initial denaturation at 94 °C, 10 cycles of 
30 s at 94 °C, 45 s at 58 °C decreasing by 0.5 °C each cycle, 2 
min at 72 °C, and 22 cycles of 30 s at 94 °C, 45 s at 52 °C, 

2 min at 72 °C. After the PCR, the products were purified using 
modified Exo-SAP method [16] with 10 U of Exonuclease I 
(Fermentas) and 0.5 U of FastAP (Fermenats). The PCR 
products were sequenced with the primers used in 
amplification reactions and with the following internal primers - 
cytB: CytB-lntF [17] and CytB-lnt2R (5'- 
ATTTGACCCTGTTTCGTGGA-3'), ND2: ND2-lntF (5'- 
ATTCAAACAGCACAAACCAT-3') and ND2-lntR (5'- 
TCGTAGTTGGGTTTGATTTA-3'). Sequencing reactions were 
performed using BigDye3.1 chemistry (Applera). Sequencing 
products were sequenced with AbiPrism 3130x1 Genetic 
Analyser (Applera). Chromatograms were preanalyzed with 
Sequencing Analysis Software ver. 5.3 and assembled in 
contigs using SeqMan program from DNAStar package 



PLOS ONE | www.plosone.org 



3 



December 2013 | Volume 8 | Issue 12 | e82464 



Phylogeography of the Carpathian Barbel 



(Lasergene). Sequences were aligned and trimmed to the 
shortest available sequence using Bioedit ver. 7.1.1 1 [18]. Due 
to technical problems not related to DNA polymorphism, we 
only obtained sequences from the 544 bp fragment of ND2 
gene. For the ATP locus we obtained a fragment of 842 bp in 
all sequences while in the cytB gene the entire amplified 
fragment of 1193 bp was successfully sequenced. All the 
unique haplotypes have been submitted to the NCBI GenBank 
nucleotide database under following accession numbers: 
KF819397-KF819405 and KF826494-KF826507. 

Previous studies demonstrated that B. carpathicus possess 
relatively low genetic variation as compared to its congeners 
[5], [6], [7]. Only seven primer pairs out of forty-five tested 
amplified polymorphic microsatellite loci in the Carpathian 
barbel: Barb37, Barb59, Barb79 [19], MFW17, MFW19, [20], 
Lco4 [21], LC293 [22] (Table 1). Loci were amplified in 
MJResearch PTC-200 Thermal Cycler in 10 pi reactions 
containing 1x HotStar Master Mix Kit (Qiagen), 200 nM of each 
primer, and 20 ng of template DNA, using the following thermal 
profile: denaturation/enzyme activation - 15 min at 95 °C; 32 
cycles of: 30 s at 94°C, 40 s at the locus optimal annealing 
temperature (Table 1), 50 s at 72 °C; final extension - 5 min at 
72°C. The efficacy of amplification in varying PCR conditions 
was validated on 3 DNAs of B. carpathicus from the Somes, 
Bela Orava and Dunajec rivers. Product lengths were 
estimated in GenMapper 4.0 (Applera) based on 
chromatograms from AbiPrism 3130x1 Genetic Analyser 
(Applera) by comparison to Rox-400HD size standard 
(Applera). 

Statistical analyses 

Genetic variation within populations. Levels of genetic 
variation were estimated overall and within populations using 
both mitochondrial sequences and microsatellite allele 
frequencies. Because the sequence variation in mitochondrial 
genome was low, all three regions were concatenated. For the 
2,579 bp fragment of mtDNA basic statistics such as number of 
haplotypes (h), number of segregating sites (S), nucleotide and 
haplotype diversities (tt and H d ) were calculated in DnaSP ver. 
5.10.01 [23]. The evolutionary relationships between 
haplotypes were estimated by maximum parsimony using TCS 
ver. 1.21 [24]. 

The expected heterozygosity (H E ) and Shannon index of 
diversity (SI) were calculated from microsatellite data with MSA 
ver. 4.05 [25], whereas allelic richness (R s ) and inbreeding 
coefficient (F, s ) were calculated in Fstat ver. 2.9.3.2 [26]. The 
Ipel river population, represented by only 2 individuals, was not 
included in R s calculation. Mean values across all loci within 
each population were calculated. Because the location of the 
amplified loci is unknown we have tested the genotype 
frequencies for signs of linkage disequilibrium using Fstat. In 
each population we performed locus by locus exact test of 
Hardy-Weinberg equilibrium as implemented in Arlequin ver. 
3.5.1.3 [27] with 1,000,000 steps in Markov chain and 100,000 
dememorization steps. 

Demographic history of the species. To check for a sign 
of recent population expansion in the variation of mtDNA 
sequences, Tajima's D and Fu's F s were calculated in Arlequin 



and tested for significance using 10,000 simulations. The 
pattern of expansions was further tested based on the 
distribution of pairwise differences between haplotypes [28]. 
Mismatch distributions were analyzed in Arlequin using both 
demographic and spatial expansion models [29]. The statistical 
support of the two models was tested using the sum of squared 
deviations (SSD) with 10,000 bootstrap replicates. 

Time to the most recent common ancestor (tMRCA) of the 
mitochondrial genome was estimated by means of Bayesian 
inference as implemented in Beast ver. 1.7.5 [30]. The 
Bayesian Skyline coalescent tree prior was selected as the 
most likely demographic scenario for the species - results of 
other tests (Tajima's D, Fu's F s , mismatch distribution, see the 
Results) suggest that the species has undergone expansion 
since the last glacial period, so constant population size in this 
period is very unlikely. On the other hand, clear limits for further 
expansion occur in the study area, thus simple expansion tree 
priors would poorly reflect the actual demographic history of the 
species. Because there are no long branches in the haplotype 
tree a simple skyline model with constant population size 
between coalescent points was used. First, short runs of fifty 
million generations sampled every 1000 generations were 
performed to choose the most likely molecular clock model and 
sequence substitution models. The resulting log files were 
inspected in Tracer ver. 1.5 (Rambaut and Drummond, http:// 
tree.bio.ed.ac.uk/software/tracer/ ). Twenty million generations 
at the beginning of each run were trimmed as the burn-in 
period after which chains reached stationarity. Based on the 
analysis of Bayesian factors, the exponential relaxed clock 
model was selected while no important differences were 
detected between different substitution schemes, thus the HKY 
model without data partitioning and site heterogeneity was 
selected as the simplest model available in the program. After 
completing preliminary runs, 5 runs of 500,000,000 generations 
sampled every 10 000 generations were performed using the 
selected model. The chains, after trimming the burn-in period 
required to reach stationarity, were combined in LogCombiner, 
ver. 1.7.5. Time estimates were obtained using the mean (0.07 
changes/site/Myr) and standard deviation (0.032) calculated 
from mutation rate estimates from recent geological events 
(isolation events 6-9 from Table 1 in [31]). A Bayesian Skyline 
Plot was constructed to show the demographic history of the 
species and the tree root height was used as an approximation 
of tMRCA of the whole mitochondrial variation. 

Genetic structure. For all pairs of populations F ST [32] was 
calculated in MSA ver. 4.05. Statistical support was estimated 
based on 10,000 permutations. Factorial correspondence 
analysis (FCA) was performed in Genetix ver 4.05.2 [33]. 

An assignment test was performed in Structure ver. 2.3.4 
[34]. Large pairwise F ST values between populations from 
distinct river basins (see the Results) suggest that current 
levels of gene flow are relatively small, thus, simpler, less 
parameterised model without admixture was selected. The 
model is appropriate to represent populations between which 
rates of gene flow are low and is more powerful at detecting 
subtle structure [34]. Initial runs of the program showed that 
clustering is very poor due to the low level of variation in the 
analyzed loci, therefore, LOCPRIOR was used to assist 
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clustering [35]. Correlated allele frequencies between 
populations were assumed to minimize spurious clustering 
between populations. Markov Chain Monte Carlo simulations 
were run for 10 6 generations with an initial burn-in period of 10 5 
generations. Runs were repeated 10 times for K ranging from 1 
to 10. The method of Evanno [36] was used to detect the most 
probable number of clusters as implemented in the 
StructureHarvester website [37]. For selected simulations the 
cluster membership probabilities were estimated in Clumpp 
[38]. 

Although, Hubisz et al [35] demonstrated that using 
LOCPRIOR does not lead to discovery of false population 
structuring, robustness of the assignment was confirmed in 
TESS ver. 2.3.1 [39]. The Bayesian clustering algorithm 
implemented in TESS incorporates geographical location as an 
additional prior [40]. The analyses were repeated 10 times for 
K=2 and 3, using model without admixture. The algorithm was 
run for 100,000 sweeps with 10,000 sweeps as the burn-in 
period. 

Results 

Mitochondrial DNA 

Within the sequenced 2,579 bp of mitochondrial DNA (-15 % 
of a complete mitochondrial genome [41]) analysed in 78 
individuals we found 19 haplotypes (Figure 1) and 22 
polymorphic sites. The largest pairwise difference found 
between haplotypes was 5 mutations (0.2 % divergence). The 
mean time to tMRCA of all the haplotypes of B. carpathicus 
was estimated by Beast to correspond to 7.8 thousand years 
BP (95% HPD: 19.3-1.15 kya). The populations from the three 
major basins were not reciprocally monophyletic (Figure 1). 
The overall most frequent haplotype was H1, and although it 
was found in all 3 river basins, H1 (and its mutational variants) 
clearly predominated in the river basins of the Vistula 
tributaries. The second most abundant haplotype H2 was most 
frequent in the Danube basin, except in the tributaries of the 
Upper Tisza river (the Topl'a, Laborec, Uzh, Iza and Somes), 
where the populations carried H3 and its derivates. 

Microsatellite loci 

Out of 45 primer pairs tested, only 7 pairs amplified 
polymorphic loci in B. carpathicus (Table 1). Primers MFW17 
and Lco4 amplified duplicated (encoded as "a" and "b"), 
polymorphic loci scorable in B. carpathicus, while primers 
LC293 probably amplified duplicated loci with only one locus 
polymorphic - the allele in the monomorphic locus possess a 
similar stutter band pattern as the alleles in the polymorphic 
locus. The number of alleles found in B. carpathicus varied 
between 3 (Lco4b) and 41 alleles (Barb59) per locus. Allelic 
richness (R s ) and expected heterozygosity (H E ) were lowest in 
Lco4a (1.6 and 0.13 respectively) and highest in Barb59 (6.4 
and 0.87 respectively). Genotype frequencies neither show 
signs of linkage disequilibrium nor deviate significantly from 
Hardy-Weinberg expectations, except for the population 
sample from the Hron river basin (including sample from the 
Sikenica river) that showed strong heterozygote deficiency (F, s 
= 0.406, Table S1). 



Genetic variation within populations 

There were substantial differences in genetic variation levels 
among the populations studied (Table 2, Figure 2). Both in the 
mitochondrial sequences and in the microsatellite allelic 
variation the populations from the Danube basin were more 
variable (h = 14, S = 17, tt = 7.3 *10 3 , H E = 0.48, R s = 3.02) 
than those from the Vistula and its tributaries (h = 8, S = 8, tt = 
3.1 *10 3 , H E = 0.40, R s = 2.55). The highest number of private 
haplotypes (8) was found in the Tisa basin, while in the Vistula 
there were 5 private haplotypes. No private haplotypes were 
found in the Strwiaz stream (the Dniester basin), however, only 
2 haplotypes from that population were obtained. 

The highest level of variation among all populations was 
found in south-eastern part of the species range (the Somes 
and Iza rivers), while the lowest values characterised the 
populations from the Strwiaz (the Dniester river basin). A very 
low level of variation was also observed in the San, Mala Wisla, 
and Orava basins, dines of variation could be observed across 
the tributaries of the Tisa/Danube with the levels of variation 
decreasing westwards, and in the Vistula river system with the 
variation decreasing both east- and westwards and with the 
populations from the Dunajec basin being most variable. 

Genetic structure 

Despite the low levels of variation, the species shows 
surprisingly strong genetic structuring. In 351 pairwise F ST 
estimates only 39 values (11.1%) did not differ significantly 
from 0, while 125 (35.6 %) were higher than F ST = 0.2 (Table 
S1). The highest value was found between the populations in 
the Strwiaz and Ipel (F ST = 0.717) and the lowest was between 
populations in the tososina and Dunajec (-0.005). In general, 
F ST values among populations from the Vistula, Danube, and 
Dniester river basins were higher than among population within 
those basins, but with several notable exceptions. The 
estimated F ST between the Orava drainage and remaining 
populations within the Danube basin were higher (mean F ST = 
0.311 ±0.125) than with populations from the other two basins 
(mean F ST = 0.157±0.0.56), and the lowest F ST among the 
Orava basin and other populations was found in comparison 
with the Mata Wisla, Sola, Skawa, Skawinka, and with the 
Wisloka rivers. The lowest F ST values found between rivers of 
the Vistula basin and rivers from other basins were between 
the pairs Poprad-Hornad (F ST = 0.051), tososina-Uh (F ST = 
0.055), and tososina-Somes (F ST = 0.063), and these values 
were lower than in comparison with other rivers from the 
Vistula basin other than those of the Dunajec river system and 
the Uszwica river. The same pattern can be observed in 
factorial correspondence analysis (Fig. 3). The first axis 
explaining 27.28 % of variation, places south-eastern 
populations and populations from the Dunajec, tososina, 
Poprad and Uszwica rivers in the centre, remaining Vistulan 
populations and the Orava basin on one side and populations 
from the south-western part of the range on the opposite side. 
The second axis groups south-eastern populations at the top of 
graph and the remaining ones at the bottom. 

An analysis of K based on runs of Structure indicated that 
the most likely number of clusters is K = 2 (K = 314), but the 
results for K = 3 and K = 6 are also presented as their K was 
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Table 2. Levels of genetic variation in 2,579 bp region of mitochondrial DNA and nine nuclear microsatellites across 27 
populations of Barbus carpathicus. 



No. 


Population 


Nmt 


Hd 


7T 


Nms 


H E ± S.D. 


R S ±S.D. 


SI ± S.D. 


1 


Somes 


2 


1 


0.00155 


11 


0.62 ±0.186 


3.83 ± 1.96 


1.18±0.60 


2 


Iza 


1 


- 


- 


8 


0.64 ±0.214 


3.87 ± 2.48 


1.26 ± 0.76 


3 


Uh 


5 


1 


0.00116 


27 


0.57 ± 0.250 


3.45 ± 2.21 


1.15±0.79 


4 


Laborec 


2 


1 


0.00039 


36 


0.56 ± 0.245 


3.24 ± 1.98 


1.12±0.75 


5 


Topl'a 


2 


1 


0.00039 


25 


0.51 ± 0.244 


2.90 ± 1.90 


0.95 ± 0.69 


6 


Torysa 


- 


- 


- 


42 


0.51 ±0.281 


3.31 ± 2.22 


1.08 ±0.85 


7 


Hornad 


6 


0.6 


0.00034 


18 


0.56 ±0.196 


3.24 ± 1.70 


1 .04 ± 0.60 


8 


Slana 


3 


1 


0.00078 


14 


0.48 ± 0.258 


3.08 ± 1.91 


0.90 ±0.67 


9 


Morgo 


2 


0 


0 


5 


0.37 ± 0.361 


2.50 ±2.00 


0.59 ± 0.63 


10 


Ipel 


1 


- 


- 


2 


0.35 ± 0.348 


N/A 


0.39 ± 0.40 




Hron 


3 


0.67 


0.00052 


10 


0.33 ± 0.299 


2.01 ± 1.27 


0.54 ± 0.54 


12 


Orava drainage 


6 


0.33 


0.00026 


18 


0.30 ± 0.233 


1.81 ±0.76 


0.46 ± 0.37 




Danube basin 


33 


0.84 


0.00073 


216 


0.48 ±0.117 


3.02 ± 0.67 


0.89 ± 0.31 


13 


Mala Wista 


3 


0 


0 


30 


0.31 ± 0.267 


2.13 ±0.89 


0.52 ± 0.47 


14 


Sota 


5 


0.4 


0.00016 


48 


0.33 ± 0.293 


2.20 ± 1.16 


0.59 ±0.57 


15 


Skawa 


1 


- 


- 


24 


0.37 ± 0.262 


2.36 ± 1.17 


0.66 ±0.51 


16 


Skawinka 


5 


0.4 


0.00031 


28 


0.37 ± 0.256 


2.33 ± 1.12 


0.67 ± 0.50 


17 


Raba 


3 


0 


0 


30 


0.38 ±0.296 


2.67 ± 1.65 


0.75 ± 0.67 


18 


Uszwica 


1 






12 


0.46 ± 0.284 


2.87 ± 1.89 


0.84 ± 0.65 


19 


Lososina 








31 


0.55 ± 0.206 


3.23 ± 1.75 


1.07 ± 0.65 


20 


Dunajec 


3 


0.67 


0.00078 


19 


0.49 ± 0.231 


3.01 ± 1.78 


0.92 ± 0.65 


21 


Poprad 


2 


1 


0.00078 


54 


0.54 ± 0.235 


3.33 ± 1.99 


1.11 ± 0.78 


22 


Biata 


3 


1 


0.00103 


68 


0.46 ± 0.234 


2.81 ± 1.58 


0.89 ± 0.64 


23 


Wistoka 


6 


0.6 


0.00026 


76 


0.46 ± 0.241 


2.62 ± 1.41 


0.84 ± 0.57 


24 


Wistok 


2 


0 


0 


42 


0.28 ± 0.294 


2.10± 1.30 


0.51 ±0.58 


25 


Lower San 


9 


0.22 


0.00009 


61 


0.30 ± 0.298 


2.04 ± 1.16 


0.53 ± 0.57 


26 


Upper San 








70 


0.30 ± 0.281 


1.99 ± 1.07 


0.52 ± 0.52 




Vistula basin 


43 


0.41 


0.00031 


593 


0.40 ± 0.093 


2.55 ± 0.45 


0.74 ± 0.20 


27 


Strwiaz (Dniester basin) 


2 


0 


0 


32 


0.16 ±0.219 


1.73 0.96 


0.31 ± 0.43 



N m t - number of haplotypes analysed, Nms - 
doi: 10.1371/journal.pone.0082464.t002 



number of individuals typed at microsatellite loci. Further explanations in text. 



notably higher (8.0 and 5.7 respectively) than for other Ks (K < 
1). For K = 2, most samples from the Danube basin strongly 
assign to a single cluster (0 2 ), except for the Orava basin 
assigned to the second cluster (B 2 ), and several individuals 
from populations in the Topl'a and Laborec rivers with 
uncertain clustering (Fig. 4). In the Vistula basin, populations 
from the Dunajec, Poprad, Lososina, and Uszwica rivers were 
assigned to the 0 2 cluster while populations from the Mata 
Wista, Sota, Skawa, Skawinka, Raba, Wistoka, Wistok, San 
rivers along with population of the Strwiaz river from the 
Dniester basin to the B 2 cluster. The assignment of the 
individuals from the Biata river was uncertain. When K = 3 was 
assumed, a major cluster from the Danube basin (Y 3 ) was 
found in most populations except for the Orava basin, Topl'a, 
and Laborec rivers, which also contained individuals from the 
Dunajec, Poprad, and Uszwica rivers. Individuals from the 
Topl'a and Laborec rivers were assigned to the same cluster 
(0 3 ) as from the Wistoka and Biata rivers, while the third cluster 
(B 3 ) unambiguously gathers population from the Orava basin 
(the Danube basin), Mata Wista, Wistok, San and Strwiaz rivers 
(the Dniester basin). Populations from the Raba, Skawinka, 



Skawa and Sota rivers were clustered either to 0 3 or B 3 . In the 
third test with K = 6, the only cluster divided between river 
basins was cluster G 6 found in all 3 river basins (the Orava 
drainage, the Mata Wista and Strwiaz rivers). The remaining 
Danubian populations were assigned strongly to either of two 
clusters (Y 6 and V 6 ). In the Vistula basin the Dunajec, Poprad, 
Lososina and Uszwica populations formed cluster P 6 , 
populations from the Wistok and San rivers were grouped into 
cluster S 6 , all individuals from the Wistoka river were assigned 
to clade O e , and individuals from other populations were 
assigned to two (Biata - P 6 +B 6 ) or more clusters (Raba - 
G 6 +P 6 +B 6 +0 6 , Skawinka and Sota - G 6 +0 6 +B 6 , Skawa - 
G 6 +0 6 +P 6 ). 

Analyses in the TESS program were performed for the most 
likely K = 2, and K = 3. For K = 2, populations from the Danube 
river basin were clustered together except for one population in 
the Orava basin that was clustered with most populations from 
northern slopes of the Carpathians except for populations in 
the Uszwica, Lososina, Dunajec and Poprad rivers (Figure 5). 
For K = 3, an additional cluster included populations from the 
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Figure 2. Levels of genetic variation in the 27 populations of Barbus carpathicus. 

doi: 10.1371/joumal.pone.0082464.g002 



0 



Laborec, Topl'a in the Danube basin and the Wistoka, and one 
population in the Biata in the Vistula basin. 

Historical demography 

Signs of population expansion were inferred from the 
majority of the tests (Table 3, Figure 6). Tajima's D and Fu's Fs 
have statistically significant negative values both in the whole 
species and in the two major river systems analysed separately 
- the Danube and the Vistula. Expansion patterns tested with 
mismatch distribution were different depending on the 
population studied. In the whole species, the mismatch 
distribution analysis did not yield statistically significant results. 
The population in the Danube basin showed moderate signs of 
both demographic and spatial expansion (p = 0.01), while the 
population from the Vistula basin showed strong evidence for 
demographic expansion (p < 10 5 ) while analysis of spatial 
expansion was not statistically supported (p = 0.72). The 
Bayesian Skyline Plot also shows a strong expansion trend 
since approximately 4000 years BP when the whole species 
appears to have reached its lowest effective population size 
(Figure 6). The expansion trend stabilizes a few hundred years 
BP, with a slightly negative tendency at the beginning of the 
coalescent simulation. 



Despite using various model settings, all attempts to run the 
coalescent simulations for smaller subsets of sequences failed. 
Probably due to low variation, the chains generated in these 
runs mixed poorly and some of the parameters were trapped in 
local peaks of probability or optimization tended towards values 
below or above the limits of variables handled by the program. 

Discussion 

Glacial refugia and colonization routes 

All the variation presently found in the Carpathian barbel 
apparently was accumulated only since the last glacial 
maximum (LGM, -22 kya). Most variation that could have 
accumulated in the species since its divergence from other 
species was thus erased during the late Pleistocene/early 
Holocene which invokes some dramatic demographic collapse 
across all the species range. The species probably survived 
the Pleistocene in glacial refugium located in the Upper Tisza 
river system [7] which is supported by the highest levels of 
genetic variation found in populations in the northern 
Transylvania. Mean tMRCA estimated with an upper bound 
located right after the LGM would appear to suggest that the 
Holocene was a turning point in the recent history of 6. 
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Figure 3. Results of factorial correspondence analysis (FCA) performed in Genetix, showing multivariate relationships 
among populations using microsatellite variation. 

doi: 10.1371/joumal.pone.0082464.g003 



carpathicus. Due to a complete lack of variation from the earlier 
history of the species, we cannot answer if this is a result of a 
single demographic collapse event or an effect of prolonged 
genetic drift during the Weichsel glaciation (110-10 kya). 
Because long stretches of lowland rivers act as effective 
barriers in rheophilic fishes, it is likely that after the colonization 
of northern slopes of the Carpathians, B. carpathicus became 
"trapped" in this isolated, remote location, and was subjected to 
more extreme demographic fluctuations than its Central 
European, rheophilic congeners - B. balcanicus and B. petenyi 
[6], [9]. Taking into account that LGM glaciers in the 
Carpathians had very narrow altitudinal range [42] and that the 
existence of Carpathian refugia were documented in a number 
of species (e.g. [43], [44]), it is possible that appropriate 
conditions for B. carpathicus could have existed in the central 
part of the Carpathian range during the entire period of 
glaciation. On the other hand, the end of the glacial period was 
accompanied by droughts that might have altered dramatically 
the functioning of riverine ecosystems, however, our data does 
not allow us to test those hypotheses. 

The current distribution of the Carpathian barbel can only be 
explained with an assumption that the species was able to 
pass segments of lowland rivers and cross watersheds several 



times. After the LGM, it is likely that B. carpathicus began to 
expand its range and crossed the Black Sea-Baltic watershed 
in a number of stream captures. Such a mode of dispersal of 
montane fishes has been documented in other parts of Europe 
in the case of the bullhead, Cottus gobio [45], [46] [47], and 
may be regarded as effective [31], [48], [49]. The range 
expansion probably comprised of several episodes after which 
the newly established population grew in number and served 
as a source of individuals for subsequent colonization [50]. Due 
to low genetic variation, a detailed picture of expansion cannot 
be unambiguously drawn, however, major patterns are well 
supported in our results. An assignment test performed in 
Structure and TESS show a few well defined clusters. Within 
the Danube basin the populations from the south-western part 
of the range (the Torysa, Hornad, Slana, Ipel, Morgo, Hron) 
consequently assign together independent of assumed K 
number. These populations also form a separate cluster in FCA 
analysis and are most likely derived from a single ancestral 
population. Populations from the central part of the Tisza basin 
do not assign consequently. There can be several explanations 
for this pattern. Firstly, this part of the range could have been 
colonized in several expansion episodes. Secondly, if the 
population was genetically diverse, a subsequent demographic 
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Figure 4. Proportional membership (q) of all individuals 
genotyped at microsatellite loci, for 3 different number of 
clusters (K = 2, 3, 6) identified by Structure. 

doi: 10.1371/joumal.pone.0082464.g004 



bottleneck could cause diversification among descendant 
subpopulations. 

Our results indicate that expansion to the northern slopes of 
the Carpathians was accomplished in at least two or three 
independent episodes. Populations from the Mata Wista, 
Strwiaz and Orava basin are genetically related and possess 
very low genetic variation, which suggest these three 
populations are remnants of the first expansion event to the 
northern slopes of the Carpathian range. Unfortunately, neither 
sequence genealogies in a long segment of mtDNA nor private 
allele distribution could shed light on the direction of expansion 
among those populations. Their genetic uniformity may suggest 
that either the founder population was genetically uniform or 
experienced a strong bottleneck right after crossing the 
watershed. The source of the colonizing individuals and the 
direction of expansion cannot be unambiguously identified. 
Variation within mitochondrial haplotypes greatly overlaps 
between the Vistula basin and south-western Slovakian 
populations of B. carpathicus, though both distance method 
(Table S1) and assignment tests using nuclear loci (Figures 4 
and 5) show those populations as genetically distinct. The 
lowest F ST values between those three populations and the 
populations from the Danube basin were found in comparison 
with populations from the Topl'a and Laborec rivers, which may 
suggest that this region was the source of the colonizing 
individuals. The second colonization event can be inferred from 
assignment tests for K=3 in Structure and TESS and includes 
populations from the Biata and Wistok rivers, that again assign 
to populations from the Topl'a and Laborec rivers. The third 
wave of expansion took place between the Hornad and 
Dunajec basins (probably from the Hornad to the Poprad). 
Those populations are closely related while being different from 
the neighbouring populations. The F ST values between them 
are some of the lowest observed among all pairwise 
comparisons (Table S1). 

Origin of genetic dines 

When dines of genetic variation are observed along the 
possible postglacial expansion routes it is often assumed that 
they arise as a result of several consecutive bottlenecks [12], 
[51]. In our study, this pattern seems more complex. Clines of 
variation were observed both in the Danube and Vistula basins. 
In the Danube basin, indices of genetic variation decrease 
gradually westwards, while in the Vistula basin the highest 
variation was observed in the Dunajec basin and it declines 
both eastwards and westwards. The origins of the clines may, 
however, be very different. The cline of variation in the Danube 
basin (except for population from the Orava basin) can be 
attributed to expansion through watersheds or long stretches of 
lowland rivers. Barbus carpathicus is strongly rheophilic and 
longer lowland river segments form effective dispersal barriers 
[8], which is supported by high F ST values between populations 
separated by river reaches with slow current (Table S1). 
Expansion through those barriers was probably accomplished 
by a limited number of individuals. Once the population is 
established in the appropriate habitat it grows in numbers until 
approaching the carrying capacity of the new site [50], which 
scenario is particularly likely in a cluster of populations from the 
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Vistula 




Figure 5. Spatial distribution of clusters detected in an assignment test performed in TESS. 

between clusters when K = 2, dashed line - borders of the 3 rd cluster for K=3. 

doi: 10.1371/joumal.pone.0082464.g005 



Solid green line - border 



south-western part of the Carpathian barbel range, that 
probably originated from a single population as evidenced by 
an assignment test in Structure but shows declining levels of 
variation from east to west. This result shows that dines of 
genetic variation may occur not only over large geographical 
scales [12], [50] but also among local populations isolated by 
effective barriers to gene flow. In the Vistula, the scenario has 
most likely been more complex. The first colonization probably 
established large, but genetically highly uniform populations, 
which resulted in the genetic proximity of currently vicariant 
populations (the Mata Wiste, Strwiaz, Orava basins). In 
contrast, subsequent expansion from the Hornad to Poprad 
and then to all basin of the Dunajec and the Uszwica probably 
consisted of a large group of genetically more variable 
individuals. Populations from the Hornad and Dunajec basins 



are genetically related and posses a similar level of genetic 
variation. A westward cline (the Raba, Skawinka, Skawa, Sola, 
Mata Wista) shows, however, no significant admixture from the 
Dunajec and Uszwica, while some genetic admixture from 
populations from the Wistok/San and the Wistoka/Biata groups 
is visible from Structure results (Figure 4). The origin of the 
eastward cline is more difficult to explain. The population from 
the Biata show traces of admixture with individuals from other 
tributaries of the Dunajec, but with a strong genetic fraction 
shared with the Wistoka population in the results from both 
Structure and TESS programs. Individuals from the Wistok and 
San assign consequently to a single population, genetically 
related to the first wave of expansion. It is possible that this 
cline was formed by independent events and should not be 
considered as a true single cline. Such a hypothesis 
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Figure 6. Bayesian Skyline Plot of population dynamics of the Carpathian barbel. Thick curves represent mean inferred 
effective population size, dashed curves mark the 95% highest probability density (HPD) intervals. Vertical line - mean 
mitochondrial tMRCA estimate of the whole species. 
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Table 3. Results of Tajima's D, Fu's F s and mismatch 
distribution tests of Barbus carpathicus expansion along 
with their statistical support. 



Demographic 
expansion 



Spatial 
expansion 



Population 


Tajima 




Fu 




model 




model 






D 


P 


Fs 


P 


SSD 


P 


SSD 


P 


Whole 
range 


-2.05 


0.002 


-14.279 


<10" 5 


0.051 


0.114 


0.036 


0.226 


Danube 


-1.84 


0.012 


-8.31 


<10" 5 


0.065 


0.01 


0.064 


0.01 


Vistula 


-1.61 


0.036 


-4.31 


0.006 


0.224 


<10" 5 


0.001 


0.72 



doi: 10.1371/joumal.pone.0082464.t003 

contradicts the previous findings of Konopihski et al. [5], which 
explained the observed variation in a part of the Vistula basin 
by expansion associated with a loss of variation. It is possible 
that such 'admixture' dines also exist in other species and that 



only range-wide studies such as this one can elucidate the true 
origins of the observed patterns. 

Demographic changes 

The Bayesian Skyline Plot shows that the demographic 
expansion of B. carpathicus started approximately 2-6 kya, 
which falls well into the mid-Holocene (Fig. 6). Demographic 
expansion is also evidenced by other genetic tests (Table 3). It 
appears that genetic drift preceding the expansion has erased 
all the variation that accumulated in the species since it 
diverged from its congeners in the Miocene. Our results may 
suggest that after the early bottleneck, the species 
encountered other demographic fluctuations. The genetic 
variation carried by individuals from the third colonization wave 
(from the Hornad to Dunajec and Uszwica) strongly altered the 
genetic composition of populations within the Dunajec and 
Uszwica river systems. It is very likely that between the first 
and third colonization episodes, the Dunajec and its tributaries 
had held populations of the Carpathian barbel - the Dunajec 
basin is the largest river system among the Carpathian 
tributaries of the Vistula, and offers a number of sites with 
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habitats suitable for B. carpathicus [8]. The Dunajec was the 
most likely colonization route between the Vistula and the 
Orava basins - the Carpathian watershed is very low between 
the Dunajec and the Czarna Orawa and passes the flat area of 
the Podhale and northern Orava regions, that is surrounded by 
the high ranges of the Tatra and the Beskid Zywiecki 
mountains. In fact, populations from the Dunajec take 
intermediate positions on FCA plot (Fig. 3) suggesting that the 
Dunajec is admixed with genes from the Hornad. If the earlier 
waves of expansion were genetically uniform, one could 
hypothesise that the successful expansion of genes carried by 
individuals from the Hornad population could be explained by 
some selection-driven processes such as heterozygote 
advantage or positive selection. A similar process may have 
shaped variation in the Biata population, which appears to have 
retained a significant fraction of distinct genetic variation 
originating from the second colonization wave. If no selective 
pressures were involved, the Dunajec and its tributaries 
(except for the Biata) should be unoccupied or contain small 
populations that were overwhelmed by individuals from the 
Hornad. Also population from Uszwica clusters unambiguously 
with populations from the Dunajec, tososina, Poprad and at 
lower K's with populations from the south-eastern part of the 
species' range. Both FCA plot and F ST values indicate that this 
population is more genetically similar to the population from the 
Hornad, thus, it is tempting to hypothesise the Uszwica was 
recolonized immediately after the third colonization episode, 
when the colonizing population was less intermixed with the 
remnant, original populations from the Dunajec river system. 

Demographic breakdowns suggested by genetic data would 
imply that B. carpathicus has very strict ecological 
requirements and is sensitive to habitat alterations. The 
species is present only within the Carpathian range and did not 
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